rm(list = ls())
#working_folder<- "/Users/bgpopescu/Library/CloudStorage/Dropbox/covid_institutions/"
working_folder<-"//Mac/Dropbox-1/covid_paper_replication/"
setwd(working_folder)

#This is to obtain:
#-table_A3
#-figure_a4
#-figure_a5

library("readxl")
library("reshape")
library("stringi")
library("stringr")
library("dplyr")
library("tidyr")
library("stargazer")
library("corrplot")
library("Hmisc")
library("glue")
library("ggplot2")
library("glue")
library('gridExtra')
library('grid')
library('lfe')
library('broom')
library('gridExtra')
library('grid')
library('stringr')
library("plm")
library("multiwayvcov")
library("lmtest")
library("psych")
library('lfe')
library('broom')
library("broom")
library("ggrepel")
library("lemon")
library("sandwich")
library("multiwayvcov")
library("lmtest")
library("spdep")
library("nFactors")
library("FactoMineR")
library("geojsonsf")
library("xtable")

data_cov_mob = geojson_sf("./data/data_cov_mob.geojson")
data_cov_mob$pct_manufacturing<-as.numeric(data_cov_mob$pct_manufacturing)

#Step1. Turning some variables to numeric
data_cov_mob$kul<-as.numeric(data_cov_mob$kul)
data_cov_mob$nat<-as.numeric(data_cov_mob$nat)
data_cov_mob$frei<-as.numeric(data_cov_mob$frei)
data_cov_mob$soz<-as.numeric(data_cov_mob$soz)
data_cov_mob$pol<-as.numeric(data_cov_mob$pol)
data_cov_mob$inter<-as.numeric(data_cov_mob$inter)
data_nocovid_nomob<-subset(data_cov_mob, week=="00")

#Step2. Calculating the PCA index version 1
irisX<-st_drop_geometry(data_nocovid_nomob)
irisX<-subset(irisX, select = c(spo, kul, nat, frei, soz, pol, inter))
irisX<-subset(irisX, !is.na(spo))

###############################################
#Step1: Determine Number of Factors to Extract#
###############################################


ev <- eigen(cor(irisX)) # get eigenvalues
ap <- parallel(subject=nrow(irisX),var=ncol(irisX),
               rep=10,cent=.05)
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea)
#Figure A4
#Save this figure manually
plotnScree(nS)

#############################
#Step2: Tables with Loadings#
#############################

# Varimax Rotated Principal Components
# retaining 5 components 
library(psych)
fit <- principal(irisX, nfactors=3, rotate="varimax")
fit # print results

#Preparing Table_A3
x<-fit$loadings
Lambda <- unclass(x)
p <- nrow(Lambda)
factors <- ncol(Lambda)
vx <- colSums(x^2)
varex <- rbind(`SS loadings` = vx)

if (is.null(attr(x, "covariance"))) {
  varex <- rbind(varex, `Proportion Var` = vx/p)
  if (factors > 1) 
    varex <- rbind(varex, `Cumulative Var` = cumsum(vx/p))
}
bottom_table<-tibble::rownames_to_column(as.data.frame(varex), "x")
names(bottom_table)[names(bottom_table) == "x"] <- "variable"


#Preparing Table_A3
cutoff <- 0.1 # (the default for the `print.loadings()` function)
Lambda <- unclass(x)
p <- nrow(Lambda)
fx <- setNames(Lambda, NULL)
fx[abs(Lambda) < cutoff] <- NA_real_
fx <- as.data.frame(fx)
rownames(fx)
rownames(fx)[1]<-"Sport"
rownames(fx)[2]<-"Culture"
rownames(fx)[3]<-"Nature"
rownames(fx)[4]<-"Free time"
rownames(fx)[5]<-"Social"
rownames(fx)[6]<-"Political"
rownames(fx)[7]<-"Interest"
fx <- cbind(variable = rownames(fx), fx)

fx2<-rbind(fx, bottom_table)

object_latex<-xtable(fx2)

print(object_latex,
      #The following line controls where the lines are drawn
      hline.after=c(-1, 0, 7, 10),
      floating = FALSE,
      file = "./paper/tables/table_A3.tex", 
      booktabs = TRUE,
      include.rownames=FALSE)

#This to get figure A.5
result <- PCA(irisX)

